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Abstract 

A  computer  program  is  developed  to  model  a  spinning 
Intercontinental  Ballistic  Missile  (ICBM)  during  the  first 
stage  boost  phase.  The  equations  of  motion  are  derived  and 
presented  and  a  full  rotation  matrix  is  used  to  show  the 
relationship  between  a  launch-centered,  nonrotating  earth, 
inertial  reference  frame  and  the  missile  body  reference  frame. 
The  moments  of  inertia  and  aerodynamic  forces  are  derived  and 
presented.  A  feedback  controller  is  derived  which  proved  to 
be  a  necessary  addition  to  the  system  in  order  to  reduce  the 
angle  of  attack.  The  angle  of  attack  of  the  missile  produced 
adverse  effects  on  the  burnout  vector  without  the  feedback 
controller  but  the  effects  are  reduced  considerably  with  the 
controller  included.  Problem  areas  include  possibly  excessive 
nozzle  gimbal  rates  caused  by  the  feedback  controller  and  the 
need  to  change  the  initial  kick  angle  if  the  missile  is 
spinning  in  order  to  achieve  the  same  burnout  conditions  as  a 
nonspinning  missile. 


vi 


STABILITY  OF  SPINNING  ICBM  IN  FIRST  STAGE  BOOST  PHASE 


Chapter  One  Introduction 

Statement  of  the  Problem 

In  the  near  future,  Intercontinental  Ballistic  Missiles 
(ICBM)  may  be  vulnerable  to  airborne  and  space-based  lasers. 
The  most  vulnerable  stage  of  the  ICBM's  flight  is  during  the 
Stage  1  burn  when  it  is  moving  slowly  away  from  a  well  known 
position  at  the  launch  site.  If  the  ICBM  is  to  survive  this 
stage  of  flight  defensive  measures  will  have  to  be  taken.  One 
of  the  ways  to  accomplish  this  would  be  to  spin  the  missile 
about  its  symmetry  axis  to  reduce  the  dwell  time  of  a  beam  on 
a  particular  point  on  the  missile.  But,  when  this  spin  is 
introduced,  other  problems  may  arise.  The  missile's  flight 
will  describe  a  coning  motion  about  its  symmetry  axis  that  may 
cause  an  excessive  angle  of  attack.  If  the  anglv  of  attack 
remains  small  it  is  possible  that  the  position  and  velocity 
errors  at  the  end  of  the  boost  phase  may  be  unacceptably 
large.  That  is  the  problem  that  this  research  effort  will 
address . 

Introduction 

Since  the  first  days  of  computers,  modelling  of  missiles 
has  been  attempted  at  varying  levels  of  complexity.  A  review 
of  the  technical  literature  available  shows  that  everything 
from  spinning  mortar  shells  to  sophisticated  ICBMs  have  been 


modelled.  Even  though  other  specific  cases  have  been  modelled 
no  one  seems  to  have  analyzed  the  case  of  a  spinning  ICBM. 

The  remainder  of  this  chapter  presents  a  summary  of  related 
efforts . 

Synopsis  of  Related  Efforts 

In  1976  the  Aerophysics  Laboratory  conducted  research  on 
statically  stable  missiles  that  were  given  a  nominal  roll  rate 
to  average  out  lifting  effects  of  configuration  assymmetr ies 
(1).  They  showed  that  even  with  a  steady  roll,  lift 
variations  can  cause  dispersion  and  lift  nonaveraging.  This 
study  was  done  for  artillery  rounds  spinning  at  about  40 
rad/sec  after  thrust  termination.  Another  study  of  small, 
spinning  missiles  presents  the  effects  of  variable  spin  rate 
and  thrust  misalignment  on  the  natural  frequencies  and  mode 
shapes  of  the  trajectory  (2).  The  effects  were  significant 
for  small  missiles  but  there  is  some  question  as  to  how  this 
applies  to  ICBMs. 

The  second  category  of  related  work  presents  trajectory 
simulation  programs  for  space  vehicle  launchers  and  ICBMs. 

They  all  develop  the  equations  used  and  present  specific  uses 
for  their  programs  but  never  consider  the  effects  of  steady 
spin  on  the  trajectory  (3,4,5).  These  research  works,  though, 
were  good  sources  for  programming  techniques  and  missile 
property  computations. 

The  most  beneficial  reference  materials  are  the  textbooks 


on  rocket  propulsion  and  spaceflight  dynamics.  The  material 
covered  is  better  explained  than  in  research  papers  and  some 
textbooks  cover  the  case  of  spinning  missiles.  There  are 
several  good  reference  texts  available  (6,7),  but  the  most 
complete  presentation  is  given  by  Cornelisse  and  others  in 


their  book  entitled  "Rocket  Propulsion  and  Spaceflight 


Dynamics"  (8) .  This  book  is  the  main  source  of  reference  for 


this  thesis. 


Chapter  Two  Derivation  of  Equations  of  Motion 


Introduction  to  the  Derivations 


In  order  to  describe  the  motion  of  a  missile  analyti¬ 
cally,  the  equations  of  motion  must  be  derived  as  a  function 
of  time  and  numerically  integrated.  The  equations  of  motion 
can  be  set  up  as  a  set  of  ordinary  differential  equations  and 
integrated  using  a  predictor-corrector  algorithm  to  arrive  at 
a  solution  at  each  time  step.  The  predictor-corrector  algo¬ 
rithm  used  in  this  research  is  of  the  fourth  order,  meaning 
that  it  carries  along  the  last  four  values  of  the  state  vector 
and  extrapolates  these  values  to  obtain  the  next  value.  It 
then  corrects  the  extrapolated  value  to  find  a  new  value  for 

the  state  vector.  It  was  first  used  by  Haming  and  bears  his 

name  ( 9) . 

Once  the  integration  method  has  been  chosen  the  equations 
of  motion  must  be  set  up  so  they  are  compatible  with  the  in¬ 
tegrator.  One  of  the  problems  to  be  solved  is  choice  of 

reference  frames.  Some  of  the  forces  and  moments  are  more 
easily  expressed  in  a  body  reference  frame  and  others  in  an 
inertial  reference  frame.  The  solution  to  this  problem  is  to 
introduce  a  rotation  matrix  which  will  relate  the  reference 
frames  chosen.  The  two  reference  frames  are  shown  in  Figure 
2.1.  The  inertial  reference  frame,  E  =  (X,Y,Z) ,  has  the 
launch  point  as  its  origin  and  the  vehicle  reference  frame, 

£r  =  (x,y,z),  has  the  vehicle  center  of  mass,  (COM),  as  its 


ft-k-V J!Lm. ^  ^  m  **-  *  •  '*— *  ■»  *  ■»-  »  ■*-  »  •*- 
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M(dV__/dt) 


T  +  W  +  F, 


(1) 


where 


M  =  instantaneous  missile  mass 
Vcm  =  velocity  of  the  COM  of  the  missile 
T  =  thrust  of  stage  1 
W  =  gravitational  force 
Fa  =  aerodynamic  forces 
The  moment  equation  is  (8:78): 


d(I*Jl)/dt  =  -mr_  X  ( Jt  X  r_)  +  rp  X  F  +  M, 


where 

I,  =  moment  of  inertia  about  the  principal  axes 
Jl  =  angular  velocity  of  vehicle  about  the  body  axes 
m  =  mass  flow  rate 

re  =  position  vector  from  the  COM  to  the  center  of  mass 
flow 

Ma  =  aerodynamic  moments 

Eqs  (1)  and  (2)  are  vector  equations  of  motion  that  can 
be  separated  into  six  scalar  equations.  The  vectors  occurrir 
in  Eqs  (1)  and  (2)  are  resolved  like  this  (8:79): 

~cm  =  (U»V»W)E  =  velocity  of  the  COM 

A  =  (p,q,r)Er  =  rotational  velocity  about  the  COM 

T  =  (Tx,Ty,Tz)E  =  thrust  of  stage  1 

£a  =  (xa'Ya'za^£  =  aerodynamic  forces 

Ma  =  (L,  M,  N)Er  =  aerodynamic  moments 

£e  =  <xe'ye'ze)§r  =  position  vector  from  the  COM  to 
the  center  of  mass  flow 

=  (9X'%'9Z)E  =  gtavitat ional  acceleration 
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I  =  ( I  vv  ,  I,,., ,  I  _„ )  E_  =  the  moments  of  inertia 
xx  yy  zz'<«*r 

It  should  be  noted  that  second-order  terms  involving  ye,  zg, 
Ty ,  and  T z  are  neglected  because  they  are  small  compared  to 
xe  and  Tx,  respectively.  Substituting  into  Eqs  (1)  and  (2) 


yields  (8:79) : 

Mdu/dt  =  Tx  +  Mgx  +  Xa  (3) 

Mdv/dt  =  Ty  +  Mgy  +  Ya  (4) 

Mdw/dt  =  Tz  +  Mgz  +  Za  (5) 

Ixxdp/dt  =  -pdlxx/dt  +  rq(Iyy-Izz)  +  mxe(yeq+zer) 

+  L  -  mpi2p  (6) 

lyydq/dt  =  -qdlyy/dt  +  pr(Izz-Ixx)  -  mqxe2  -  xgFz 

+  zeFx  +  M  (7) 

Izzdr>/dt  =  _rdlzz//dt  +  pqUxx-Iyy)  ~  "irxe2  +  xeFy 

-  y0Fx  +  N  (8) 


where  =  offset  distance  of  the  COM  flow  from  the  missile's 
centerline 

Eqs  (3),  (4),  and  (5)  differ  from  Cornelisse's  equations 
because  he  presents  all  of  the  vectors  in  the  body  frame  and 
thus  includes  an  additional  term  necessary  to  decompose  the 
Vcm  into  the  body  frame.  That  was  avoided  in  this  treatment 
by  leaving  all  of  the  terms  in  Eqs  (3),  (4),  and  (5)  in  the 
inertial  reference  frame. 

The  last  term  in  Eq  (6)  is  not  in  Cornelisse's  equation 
but  it  is  included  here  because  otherwise  the  spin  of  the 
missile  is  not  properly  modelled  (7:225).  As  the  propellant 


is  burned  and  exhausted,  it  takes  with  it  a  small  amount  of 
ingular  momentum.  If  this  term  is  not  included  in  Eq  (6)  the 
computer  simulation  would  show  that  the  missile  spins  faster 
and  faster  as  the  propellant  is  burned  when  this  isn't 
necessarily  true. 

The  Kinematical  Equations  of  Motion 

In  the  dynamical  equations,  the  forces  and  the  moments 
are  dependent  on  the  position  and  orientation  of  the  missile. 
The  kinematical  equations  are  needed  to  relate  the  position 
and  orientation  of  the  missile  to  the  translational  and 
rotational  velocity.  The  first  three  kinematical  equations 
are  represented  by  the  vector  equation  (8:81): 


dficn/dt 

Bern  '  - 

Or,  in  component  form: 

*  ~cm 

the  position  of  the  COM 

(9) 

dX/dt 

*  u 

(10) 

dY/dt 

=  V 

(ID 

dZ  /dt 

=  w 

(12) 

Once  again,  by  leaving  the  Vcm  vector  in  the  inertial  frame, 
a  coordinate  transformation  is  avoided  that  is  included  in 
Cornelisse's  text. 

The  last  set  of  equations  necessary  is  the  transformation 
matrix  that  converts  vectors  from  the  inertial  frame  to  the 


body  frame.  They  are,  in  vector  form  (8:81) 


dAr/dt  = 


0  r  -q 
r  0  p  Ar 
q  -p  0 


where  £r  is  a  3X3  matrix  of  the  direction  cosines  between  the 
body  reference  frame  and  the  inertial  reference  frame.  Eq 
(13)  can  be  expanded  into  the  nine  equations  that  it 
represents : 


dArii/dt 

= 

rcos (y , X) 

-  qcos(z,X) 

(14) 

dAri2/dt 

= 

rcos (y  ,Y) 

-  qcos(z,Y) 

(15) 

dAr  3/ dt 

= 

rcos  (y  ,Z ) 

-  qcos(z,Z) 

(16) 

dAr  21/dt 

= 

-rcos  (x ,  X) 

+  pcos (z , X) 

(17) 

dAr  22/dt 

= 

-rcos (x  ,Y) 

+  pcos (z  ,Y) 

(18) 

dAr23/dt 

= 

-rcos (x ,Z ) 

+  pcos(z,Z) 

(19) 

dAr31/dt 

= 

qcos (x ,X) 

-  pcos (y ,  X) 

(20) 

dAr  32/dt 

= 

qcos (x ,Y) 

-  pcos (y ,Y) 

(21) 

d  Ar  3  3/^ 

= 

qcos (x ,Z ) 

-  pcos  (y  ,Z  ) 

(22) 

where  cos(i,j)  indicates  the  cosine  of  the  angle  between  i  and 
j  which  are  also  elements  of  &r 


r  =  (xfy,z)J§r  =  the  axes  of  the  body  frame 
Scm  =  =  the  axes  of  the  inertial  frame 

The  use  of  a  full  rotation  matrix  avoids  the  singularities 
associated  with  any  particular  set  of  orientation  angles.  The 
set  of  18  equations  of  motion  that  fully  describe  the  behavior 
of  the  missile  are  Eqs  ( 3 ) - ( 8 )  ,  ( 10) - ( 12)  r  and  ( 1 4 ) - ( 2 2 )  . 
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Chapter  Three  Derivation  of  Missile  Parameters  and 
Aerodynamic  Forces 

Introduction  to  Missile  Data 

All  of  the  forces  and  moments  in  Eqs  (3)  -  (8)  vary  with 
time,  position,  or  orientation  of  the  missile.  In  order  to 
find  relationships  for  each  of  these  terms  several  assumptions 
were  made.  Some  of  the  higher  order  terms  were  assumed  to  be 
small  and  thus  ignored  leaving  linear  relationships.  The  mis¬ 
sile  data  used  in  all  analyses  in  this  report  resemble  that  of 
a  Minuteman  III  (MM  III)  ICBM  but  this  is  only  for  convenience 
since  the  data  was  available  (10) .  Other  approximations  were 
made  that  will  be  presented  throughout  this  chapter. 

To  avoid  the  added  complication  of  launching  from  a  vert¬ 
ical  position  and  at  some  later  time  initiating  a  pitchover 
maneuver  for  a  gravity  turn,  a  cold  launch  system  was  assumed. 
This  is  a  departure  from  current  hardware  configurations  since 
MM  III  missiles  can't  be  cold  launched.  The  assumptions 
associated  with  the  cold  launch  are  that  the  stage  1  engine 
starts  at  50  feet  above  mean  earth  radius  with  the  center  of 
mass  moving  at  50  ft/sec  vertically.  It  is  also  assumed  the 
desired  initial  angular  velocities  have  been  achieved  when  the 
thrust  begins.  To  simulate  a  gravity  turn  trajectory  an  ini¬ 
tial  kick  angle  was  input  by  misaligning  the  body  axes  with 
the  inertial  axes  at  the  beginning  of  the  simulation. 


The  most  important  external  forces  acting  on  the  missile 


are  the  gravitational  forces  and  the  thrust.  The  gravita¬ 
tional  forces  are  dependent  on  the  distance  from  the  center  of 
the  Earth  and  can  be  represented  by  (8:75): 


/vV 


9  *  -=M£cm/Rcn,3 


where 

g  =  (9v'9u»9*)e  *  gravitational  acceleration  vector 

GM  =  1 . 407646882E+16  ft^/sec^  =  the  gravitational 
parameter  of  the  Earth 

Rcm  =  the  magnitude  of  the  position  vector  from  the 
Earth's  center  to  the  COM  of  the  missile 

£cm  =  (X»Y*Z)£  =  position  vector  from  the  Earth’s 
center  to  the  COM  of  the  missile 

This  vector  can  be  resolved  into  its  three  components,  gx, 

g^,  and  gz,  and  multiplied  by  the  instantaneous  missile 

mass,  M,  to  produce  the  respective  terms  in  Eqs  (3)  -  (5). 

The  thrust  is  assumed  to  be  a  constant  throughout  the 
duration  of  the  first  stage  burn.  This  is  never  achievable  in 
actual  engine  firings  but  the  fraction  of  a  second  that  it 
takes  for  the  thrust  curve  to  level  off  is  not  worth  including 
in  this  study.  In  order  to  use  the  thrust  in  Eqs  (3)  -  (5) 
though,  it  must  be  transformed  from  the  body  frame  to  the 
inertial  frame.  This  is  done  by  multiplying  the  transpose  of 
Ar  by  the  body  frame  thrust,  F: 


2 


where 


9 

b 

v; 

n 


i 


^  =  (Tx,Ty,Tz)E  =  thrust  in  the  inertial  frame 

Ar  =  transformation  matrix  from  inertial  to  body  frame 
F  =  (Fx,Fy,Fz)Er  =  thrust  in  the  body  frame 
T  is  then  resolved  into  its  t  ree  components  for  use  in  Eqs 
(3)  -  (5)  . 

Moments  of  Inertia  and  Their  Derivatives 

The  moments  of  inertia  and  their  time  derivatives  are 
very  important  to  this  study  because  of  their  effects  on 
rotational  velocity.  The  best  source  of  accurate  data  for 
this  would  be  experimental  results  where  the  moments  have  been 
determined  but  this  data  is  not  readily  available  in  the 
unclassified  literature.  Also,  in  order  to  use  this  computer 
program  to  simulate  other  missiles,  it  would  be  better  if  the 
program  computed  an  estimate  of  the  moments  of  inertia  from  a 
few  standard  inputs  since  the  actual  moments  are  hard  to  find. 
In  order  to  calculate  the  moments  and  keep  the  calculations 
generic,  several  assumptions  were  made.  The  first  assumption 
is  that  each  stage  is  a  solid,  homogeneous  cylinder  so  that 
the  mass  distribution  is  known.  This,  of  course,  is  not  the 
case  in  a  real  missile  but  the  results  are  probably  close 
enough  for  this  study.  The  next  assumption  is  that  the  fuel 
in  the  first  stage  motor  burns  uniformly  from  the  inside  out 
to  within  one  inch  of  the  outer  radius.  In  the  time 


derivatives  of  the  moments  of  inertia,  the  burn  rate  is  then 


considered  to  be  a  constant  given  as: 


db/dt  =  b  =  bf/burnt  (25) 

where 

bf  =  radius  of  stage  1  minus  one  inch  (for  outer  shell 
thickness) 

burnt  =  burn  time  of  stage  1 

The  final  assumption  is  that  the  missile  is  an 
axisymmetr ic ,  rigid  body  with  the  y  and  z  body  axes'  moments 
of  inertia  the  same.  Their  time  derivatives  are  also 
considered  to  be  the  same. 

There  are  several  preliminary  calculations  that  are 
necessary  for  the  final  moment  of  inertia  equations  to  be 
easier  to  understand.  The  distances  are  defined  in  Figure 
3.1. 

The  equations  are: 


m^  =  iDq  mt 

where 

m^  =  instantaneous  mass  of  stage  1 
mQ  =  initial  mass  of  stage  1 
m  =  mass  flow  rate  form  stage  1 
t  =  elapsed  time  in  seconds 

mtop  =  m2  +  m3  +  m4 
where 

mtop  =  total  mass  of  stages  2,  3,  and  4 


(26) 


(27) 


Figure  3.1  Missile  Dimensions 


xbot  = 

x9  =  LI  +  L2/2 


x3  =  LI  +  L2  +  L3/2 
x4  =  LI  +  L2  +  L3  +  L4/2 


where 


xbot  =  locati°n  of  COM  of  stage  1  from  bottom  of  missile 

X2  =  location  of  COM  of  stage  2  from  bottom  of  missile 

x3  =  location  cf  COM  of  stage  3  from  bottom  of  missile 

x4  =  location  of  COM  of  stage  4  from  bottom  of  missile 


LI  =  length  of  stage  1 
L2  =  length  of  stage  2 
L3  =  length  of  stage  3 
L4  =  length  of  stage  4 


x  top  =  <x2m2  +  x3m3  +  x4m4)/mtop  (32> 

x  =  <xbotml  +  xtopmtop)/(ml  +  mtop>  <33> 

dl  =  x  -  xbot  <34> 

d 2  =  x top  ■  x  <35> 

whete 


xtop  =  location  of  COM  of  the  top  3  stages  from  bottom 
of  missile 

x  =  distance  between  COM  and  bottom  of  missile 
d^  =  moment  arm  distance  for  stage  1 
d2  =  moment  arm  distance  for  top  3  stages 

M  =  tmass  -  mt  (36) 

where 

tmass  =  total  mass  of  the  missile  at  launch 
M  =  instantaneous  missile  mass 

The  standard  equations  for  the  moments  of  inertia  of  a  solid 
homogeneous  cylinder  are: 

Ixx  =  mass*r ad ius2/2  (37) 

Iyy  =  Izz  =  mass ( 3*radius2  +  length2)/12  (38) 


The  moments  of  inertia  for  a  four-stage  missile  are: 


=  (m1(r1^  -  b  t‘ 


+  m2r  2  +  m3r3*’ 

+  m.r42)/2 


1 2  2  =  (m1(3r12  +  LI2)  +  m2(3r22  +  L22) 
+  m3(3r32  +  L32)  +  m4(3r42  +  L42) 

-  3m1b2t2)/12  +  mtopd22  +  midi2 


where 


=  stage  1  radius 
r2  =  stage  2  radius 
r3  =  stage  3  radius 
r4  =  stage  4  radius 

Ivv  =  moment  of  inertia  about  the  body  frame  x  axis 

Iyy  =  moment  of  inertia  about  the  body  frame  y  axis 

I2Z  =  moment  of  inertia  about  the  body  frame  z  axis 

The  derivatives  of  the  moments  of  inertia  are  given  by: 


d  1  xx^dt 

dlyy/dt 


(-2m0tb2  -  mr L2  +  3mb2t2)/2  (' 

dlzz/dt  =  ("mtopd2dl  +  m1d12)2m/M 

-  ( 3r  L2  +  L 1 2 ) m/ 1 2  -  b2mQt/2 


+  3b2mt2/4  -  md12 


The  moments  of  inertia  and  their  derivatives  appear  in  Eqs  (6) 
-  (8)  . 


Aerodynamic  Forces 


The  aerodynamic  forces  and  their  moments  are  extremely 
hard  to  model  in  a  concise  way  without  oversimplifying  the 
terms.  They  depend  on  the  velocity  of  the  missile  with 
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respect  to  the  air,  the  local  atmospheric  pressure,  density, 
and  temperature  as  well  as  the  missile's  shape  and  angle  of 
attack.  The  force  equation  will  be  presented  first  and  each 
of  the  terms  in  the  equation  will  be  defined  thereafter.  The 
aerodynamic  forces,  in  vector  form,  can  be  given  by: 

Fa  -  -CdfA|Vcm|Vcm/2  (43) 

where 

F,  =  ( X, ,Y ,  ,Z  , )  E  =  aerodynamic  forces  in  the 

r*  a  a  a  a  *+  , 

inertial  frame 

Cd  =  coefficient  of  drag  of  the  missile 
p  =  local  atmospheric  density 
A  =  effective  area  of  the  missile 
Vcm  =  velocity  of  the  COM  of  the  missile 
The  coefficient  of  drag  varies  with  Mach  number, 
Reynold's  number,  angle  of  attack,  and  shape  of  the  missile, 
to  name  a  few.  The  Cd  used  in  this  report  is: 

Cd  =  Cd<>  +  Cd^c  (44) 

where 

Cd  =  drag  coefficient 
Cd#  =  base  drag  coefficient 

CH  =  drag  coefficient  variance  with  angle  of  attack 
«  =  angle  of  attack 

Most  references  on  the  subject  give  good,  accurate  descrip¬ 
tions  of  the  factors  that  influence  Cd  but  few,  if  any, 
present  the  data  required  in  this  report.  Reference  11 
presents  a  graph  (11:88)  that  shows  that  the  base  drag 


0 . 3  was 


coefficient  varies  with  Mach  number.  A  value  of  = 

ao 

obtained  from  Reference  11  (11:88).  There  was  no  data 
available  for  for  this  particular  missile  shape.  There 

are,  however,  references  to  conical  bodies  in  Krasnov  (12) 
with  Cj  =  0.1  so  that  is  the  value  used  for  this  study. 

The  variation  of  the  density  of  the  atmosphere  with 

altitude  was  simulated  using  a  simple  atmospheric  model 
(6:18-19).  It  uses  the  thermal  lapse  rate  at  varying  heights 
above  the  Earth's  surface  to  extrapolate  between  values  in  a 
lookup  table. 

The  area  of  the  missile  effecting  the  aerodynamic  forces 
is  given  by  (9) : 

A  =  Afcos<*  +  Assin*  (45) 

where 

A  =  effective  area  of  the  missile 
Aj:  =  frontal  area  of  the  missile 

As  =  surface  area  of  the  missile 

ei  =  angle  of  attack 
The  frontal  area  is  given  by: 

Af  =  IT  r  x2  (46) 

where  rx  is  the  radius  of  the  largest  stage  of  the  missile 
(stage  1).  The  surface  area  is  given  by: 

Ag  =  r  x  2  +  2TT(  r  XL1  +  r2L2  +  r -jLS  +  r^L4)  (47) 

All  of  these  variables  have  been  defined  earlier  in  this 


The  angle  of  attack  is  defined  as  the  angle  between  the 
velocity  vector  in  the  body  frame  and  the  x  axis  of  the  body 
frame.  The  velocity  vector  can  be  decomposed  into  its  three 


components  which  produces  the  three  components  of  the 
aerodynamic  forces,  namely,  Xa ,  Ya,  and  Za. 

The  aerodynamic  moments  are  given  by: 


where 


=  d 

''•a  fsi 


X  F, 

~.a 


(48) 


Ma  =  (L,M,N)Er  =  aerodynamic  moments  in  the  body 
^  f  r  a  me 

d  =  the  moment  arm  which  is  given  by: 


d  =  -(x  -  xcp)t  (49) 

where 

x  =  distance  between  COM  and  bottom  of  missile 

xcp  =  distance  between  the  center  of  pressure  and  the 
^  bottom  of  missile 

The  center  of  pressure  is  usually  determined  experimentally 
but  since  that  is  not  possible  in  this  case  the  data  from 
Krasnov  (12:442)  will  suffice.  Krasnov  reports  that  for  a 
pointed  body  of  revolution  with  a  length  of  six  to  eight 
diameters,  the  center  of  pressure  is  approximately  equal  to 
52%  of  the  missile  length,  according  to  wind-tunnel  tests.  He 
chooses  the  top  of  the  missile  as  his  zero  reference  so  48%  of 
the  total  missile  length,  28.7  feet,  will  be  used  for  this 
study.  The  length  of  this  missile  is  about  11  diameters  but 


this  difference  should  be  insignificant.  Since  the  center  of 
pressure  is  behind  the  center  of  mass  the  missile  has  an 
inherent  stability  called  "arrow  stability". 

Since  the  moment  equations  are  in  the  body  frame  the 
aerodynamic  moments  must  also  be  in  the  body  frame.  The 
transformation  is  made  by  use  of  the  transformation  matrix, 

Ar ,  and  after  the  cross  product  is  performed,  Eq  (48)  expands 
to : 

L  =  0  (50) 

M  =  -(l/2)dCdfAlVcml  (A31u  +  A32v  +  A33w)  (51) 

N  =  (l/2)dCdpA|  Vcml  (A21u  +  A22v  +  A23w)  (52) 

All  of  the  variables  above  have  been  defined  earlier  in  this 
chapter  . 

Momentum  Loss  Term 

As  the  burned  propellant  is  exhausted  out  the  nozzle  it 
takes  a  small  amount  of  angular  momentum  with  it.  A  term  must 
be  included  to  account  for  the  angular  momentum  imparted  to 
the  exhausted  propellant  or  else  the  angular  momentum  will 
remain  with  the  missile  causing  the  spin  rate  to  increase. 

This  is  not  treated  in  Cornelisse’s  presentation  but  it  is 
added  to  this  analysis.  Thomson  shows  that  this  angular 
momentum  loss  can  be  accounted  for  by  (7:225): 

Mx  =  -rhpi2p  (53) 
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where 


U*A»  J«*. 


wr?  w-’rj  wir^v 


"VW^ 


5 

/ 


Pl2 


moment  in  the  body  frame  x  direction 
momentum  loss 

mass  rate  of  flow 

the  square  of  the  offset  distance  of 
and  is  estimated  to  be  approximately 


due  to 


the  COM  flow 
(l/2)r12 


This  term  is  then  included  in  Eq  (6)  as  shown.  If  this  term 
is  not  included  the  spin  rate  rises  to  unacceptable  rates  and 
the  system  is  not  properly  modelled. 


I 


Chapter  Four  Checkout  of  Equations 
Introduction  to  Checkout 

In  order  to  ensure  that  the  computer  model  of  the 
missile's  trajectory  is  accurate,  each  term  was  checked 
individually  as  it  was  added  to  the  equations  of  motion.  In 
addition,  there  are  a  number  of  tests  that  can  be  performed  to 
make  sure  the  computer  results  match  with  hand  calculations. 
The  first  tests  involve  setting  up  driver  routines  to  make 
sure  the  numerical  integrator  works  satisfactorily.  This  was 
done  and  it  showed  that  the  numerical  integrator  gave  consid¬ 
erably  better  results  when  the  inputs  were  double  precision 
rather  than  single  precision.  For  example,  if  an  input  error 
was  made  in  the  fourth  significant  digit  the  error  was  obvious 
in  the  final  results  with  single  precision  but  it  wasn't  quite 
as  pronounced  with  double  precision  accuracy.  This  is  to  be 
expected,  of  course,  since  the  algorithm  used  extrapolates  and 
corrects  thousands  of  times  in  a  single  test  case.  As  a 
result,  double  precision  variables  were  used  throughout  the 
program. 

Torque-Free  Axisymmetric  Rigid  Body 

The  first  terms  to  be  checked  out  were  the  terms  which 
make  up  Euler's  equations  for  a  torque-free,  axisymmetric, 
rigid  body.  They  are: 

dp/dt  =  (Iyy  -  Izz)qr/Ixx  (54) 


22 


liUwM*' 


dr/dt  =  (Ixx  -  Iyy)pq/Izz 


These  equations  are  readily  identifiable  as  simplifications  of 
Eqs  (6)  -  (8).  Since  Iyy  =  Izz  for  an  axisymmetric  body, 
dp/dt  =  0  or  p  =  constant.  This  leaves  a  pair  of  constant 
coefficient,  linear  differential  equations  whose  solution  is 
known  (13).  The  magnitude  of  the  angular  velocity  vector  must 
remain  constant  and  the  period  of  oscillation  can  be  calcu¬ 
lated  by  hand  and  compared  with  computer  results.  They 
compared  favorably  to  the  fourth  significant  digit.  The  angle 
between  the  angular  velocity  vector  and  the  missile's  symmetry 
axis  must  also  remain  constant  and  this  also  proved  to  be 
tr  ue . 


Gravity  and  Position-Velocity  Relations 


The  next  terms  to  be  checked  were  the  gravity  terms  and 
the  position-velocity  relationships.  The  gravity  terms  are 
those  found  in  Eqs  (3)  -  (5)  and  represent  the  first  terms 
that  effect  the  translational  velocity  of  the  missile.  Since 
the  velocity  of  the  missile  is  the  time  derivative  of  its 
position,  Eqs  (10)  -  (12)  were  included  at  the  same  time.  The 
checkout  used  for  these  terms  was  to  set  the  missile's 
velocity  equal  to  the  orbital  speed  at  the  altitude  chosen  and 
see  if  the  missile  would  orbit  the  earth  in  the  same  period 
that  hand  calculations  predicted.  This  checked  out  to  be 


-  .  .J 


true  . 
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Transformation  Matrix 

The  transformation  equations  were  the  next  set  of 
equations  to  be  installed  and  checked.  They  are  given  by  Eqs 
(14)  -  (22).  In  order  to  check  the  equations  the  missile  was 
forced  to  cone  about  its  angular  momentum  vector .  This  was 
done  by  choosing  A^j  =  cosfl  where  0  is  the  angle  between  the 
angular  momentum  vector  and  the  symmetry  axis  of  the  missile. 
The  angular  momentum  vector  in  the  body  frame  is  given  by 
(13:16)  : 

Shod  -  lxx^  +  Jyy^  +  °*  (57) 

The  ^  component  is  chosen  to  be  zero.  The  transformation 
matrix  will  be  computed  so  that: 

Sin  =[AjHbod  “  <*g  +  +  lH^g  (58> 

where 

Hin  =  the  angular  momentum  vector  in  the  inertial 

win  J 

frame 

lg ,  5g »  =  the  unit  vectors  in  the  inertial  frame 

This  represents  a  set  of  equations  that  can  be  solved 
simultaneously  with  the  standard  equations  which  relate 
elements  of  rows  and  columns  of  a  transformation  matrix, 
namely,  the  sum  of  the  squares  of  the  elements  of  any  row  or 
column  must  equal  one.  After  solving  for  the  A^j's  and  using 
them  in  the  program,  the  Agi  element  should  not  vary  with 


time  because  it  is  a  constant  in  a  coning  motion.  That  is 
what  the  computer  results  showed  so  it  checks. 


Jet  Damping  Terms 

The  jet  damping  terms  are  caused  by  the  exhaust  jet  in  a 
spinning  missile  and  have  a  damping  effect  on  the  angular 
velocities.  In  Eqs  (6)  -  (8)  they  are  the  terms, 
respectively : 

mxe(yeq  +  zer) ,  -mqxe2,  -mrxe2 

They  result  from  the  cross  product  of  the  angular  velocity 
with  the  center  of  mass  flow  offset  vector,  re.  The  higher 
order  terms  have  been  dropped  since  they  are  negligibly  small 
compared  to  the  terms  retained.  The  magnitudes  of  the  jet 
damping  terms  are  several  orders  of  magnitude  larger  than  the 
moment  of  inertia  terms  already  in  place  so  an  approximate 
solution  can  be  used  to  check  the  jet  damping  terms'  effects. 
The  approximation  is: 

Izzdr/dt  =  Mqxe2  (59) 

Jdr/r  =  -  J  (Mxe2/Izz)dt  (60) 

r(t)  =  r (t0)exp(Mxe2(t-t0)/Izz)  (61) 

where 

r(t)  =  the  angular  velocity  about  the  z  body  axis  at 
t  i  me  t 

r(tQ)  =  the  initial  angular  velocity  about  the  z  body 
axis  at  time  t  =  tQ  =  0 

By  substituting  the  values  that  the  program  used  for  Izz,  M, 
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xe,  and  r  (tQ)  the  approximate  solution  can  be  compared  to 
Q  the  program's  results  as  shown  in  Figure  4.1.  The  results 

show  that  the  moment  of  inertia  terms  are  almost  insignificant 
and  that  the  program  is  producing  the  predicted  results. 

Figure  4.2  shows  that  the  jet  damping  terms  do  reduce  the 
angular  velocities  with  time.  In  other  words,  jet  damping  has 
the  desirable  result  of  reducing  the  coning  motion. 

Thrust 


I 


I 


The  addition  of  the  thrust  terms  to  Eqs  (3)  -  (5)  is  the 
next  step.  The  thrust  is  in  the  body  axis  x  direction  and 
must  be  transformed  to  the  inertial  coordinate  system  by  using 
the  transformation  matrix,  £r .  In  the  next  chapter,  thrust 
misalignment  will  be  used  to  control  the  angle  of  attack  but 
for  now  the  thrust  will  be  in  the  body  axis  x  direction  only. 
The  thrust  terms  were  tested  by  simulating  a  vertical  flight 
and  a  gravity  turn  trajectory.  The  simulations  worked  as 
expected  with  the  flight  path  angle  remaining  at  zero  for  the 
vertical  flight  and  growing  steadily  in  the  gravity  turn 
simulation. 

Aerodynamic  Forces  and  Moments 

The  aerodynamic  forces  and  moments  are  extremely 
difficult  to  model  accurately  and  also  very  expensive  in 
computer  time  consumption.  They  are  included  in  this  report 
with  many  simplifying  assumptions  so  that  the  computational 
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time  of  the  program  could  be  kept  lower  and  to  provide  a 
certain  amount  of  completeness.  The  terms  were  all  hand 
checked  at  various  stages  of  a  flight  simulation  and  the 
atmospheric  model  was  tested  by  checking  the  density  the 
program  output  at  various  altitudes. 

Moments  of  Inertia 

Up  to  this  point,  the  moments  of  inertia  were  estimated 
by  calculating  the  moments  of  a  solid  cylinder  with  the  radius 
of  the  first  stage  and  the  length  of  the  missile  as  its 
dimensions.  With  the  arrival  of  more  extensive  data  (10),  a 
more  accurate  model  was  constructed.  Again,  the  program's 
output  was  checked  against  hand  calculations  at  various  stages 
of  a  flight  simulation  and  the  results  were  the  same.  The 
equations  for  the  moments  of  inertia  are  presented  in  Chapter 
Three  in  Eqs  (39)  and  (40).  Table  4.1  shows  the  values  that 
were  used  in  all  the  simulations  after  these  equations  were 
implemented . 

Table  4.1  Minuteman  III  Data 


Length , ft 

Radius  ,  ft 

Initial  Mass, slugs 

Stage  1 

25.3 

2.75 

1571.14 

Stage  2 

13.1 

2.15 

482 . 38 

Stage  3 

7.1 

2.15 

253.31 

Stage  4 

14.1 

2.15 

55.32 

Totals 

59.9 

— 

2362.15 

t.' 'i*.'  -',  '7.'  y.".'.".1" 


The  stage  1  thrust  used  was  202,000  lb^  and  the  stage  1 
propellant  mass  was  1424.44  slugs. 

The  derivatives  of  the  moments  of  inertia  are  given  in 
Eqs  (41)  and  (42)  and  were  checked  in  the  same  way.  Since  the 
moments  of  inertia  are  decreasing  with  time  because  of  the 
decrease  in  mass,  the  derivatives  of  the  moments  are  negative. 
Since  they  are  subtracted  from  the  other  moments  in  Eqs  (6)  - 
(8)  their  net  effect  is  positive  or,  they  decrease  the  jet 
damping  effect.  This  is  also  the  prediction  of  Cornelisse 


inure  4.1  z-axis  Rotational  Velocity  (r)  as  a  Function 
of  Time 


Motivation  for  Feedback  Controller 

After  all  the  equations  had  been  included  and  all  the 
terms  checked  individually,  the  set  of  equations  was  analyzed. 
Some  of  the  same  characteristics  that  were  mentioned  in 
Chapter  Four  must  still  be  true  for  the  complete  system. 
Namely,  the  missile  should  fly  a  gravity  turn  trajectory  and 
the  angle  of  attack  should  drop  off  quickly  and  stay  small. 
Figure  5.1  shows  the  flight  path  angle  of  the  trajectory  as  a 
function  of  time.  The  flight  path  angle  used  throughout  this 
report  is  defined  as  the  angle  between  the  inertial  vertical 
axis,  X,  and  the  velocity  vector  in  the  inertial  frame. 

Notice  that  the  flight  path  angle  rises  to  0.59  radians  or 
about  33°  and  stays  there.  This  isn't  normal.  Figure  5.2 
shows  the  angle  of  attack  plotted  as  a  function  of  time  and 
that  the  missile  flight  simulation  is  flying  the  missile  at  a 
rather  high  angle  of  attack  throughout  the  61  seconds  of 
flight.  The  angle  of  attack  used  throughout  this  report  is 
defined  as  the  angle  between  the  missile's  symmetry  axis  and 
the  velocity  vector  in  the  body  frame.  The  initial  kick  angle 
for  the  plots  was  15°  and  the  majority  of  the  flight  the 
missile  is  flying  at  an  angle  of  attack  greater  then  0.0873 
radians  or  5°.  Figures  5.1  and  5.2  were  produced  with  only 
the  thrust,  gravity,  and  moment  of  inertia  terms  included  in 
the  simulation.  The  reason  the  angle  of  attack  increases 


during  the  flight  is  because  the  velocity  vector  falls  a way  in 
a  gravity  turn  but  the  nose  of  the  missile  doesn't  follow  it 
because  it  is  coning  about  the  angular  momentum  vector.  So 
the  feedback  controller  is  implemented  to  reduce  the  angle  of 
attack . 


Derivation 


The  general  idea  behind  the  feedback  controller  is  to 
compute  the  angle  of  attack  about  the  y  and  z  body  frame  axes 
and  to  feed  back  a  restoring  torque  to  the  moment  equations 
given  in  Eqs  (7)  and  (8).  These  two  angles  of  attack  are 
related  to  the  angle  of  attack  that  is  being  reduced  through 
spherical  trigonometry.  The  restoring  torque  is  computed  by 
first  finding  the  missile's  velocity  in  the  body  frame.  The 
transformation  matrix,  Ar ,  is  multiplied  by  the  inertial 
velocity  vector  to  find  the  components  of  the  body  frame 
velocity  vector.  Next,  the  angles  of  attack  about  the  y  and  z 
body  frame  axes  are  computed  by: 


fiy  =  arctan  ( vk/v  j. ) 
=  arctan  (v  j/v^ 


where 


=  the  angle  of  attack  about  the  y  body  frame  axis 
=  the  angle  of  attack  about  the  z  body  frame  axis 
=  velocity  component  in  the  x  body  frame  direction 
=  velocity  component  in  the  y  body  frame  direction 
=  velocity  component  in  the  z  body  frame  direction 


Next,  the  y  and  z  components  of  the  position  vector  fr 
the  missile's  COM  to  the  center  of  mass  flow  are  computed. 
This  is  done  by  (14): 

Ye  =  Kr  +  C0Z  (64) 

ze  =  Kq  +  CQy  (65) 

where 

ye  =  y  axis  component  of  the  position  vector  to  the 
COM  flow 

z0  =  z  axis  component  of  the  position  vector  to  the 
COM  flow 

K  =  gainl 

C  =  gain2 

q  =  angular  velocity  about  the  y  body  frame  axis 


r  =  angular  velocity  about  the  z  body  frame  axis 
The  gains,  K  and  C,  will  be  discussed  a  little  later  in  this 
chapter.  The  three  components  of  the  body  frame  thrust  can  be 
computed  from  the  components  of  the  position  vector  to  the  COM 
flow  using: 


F 

y 

=  Tye/|re| 

(66) 

Fz 

=  Tze/|re| 

(67) 

Fx 

=  (T2  -  Fy2  -  Fz2) X/2 

(68) 

where 

Fx,  F, 

,  F z  =  the  three  components  of 
thrust 

the  body  frame 

lre* 

=  magnitude  of  the  position  vector  to  the  COM 
f  low 

T 

=  thrust,  a  constant 

The  feedback 

terms  added  to  Eqs  (7)  and  (8) 

are : 

MFEEDJ  =  FxZe 

(69) 

MFEEDK  =  -Fxye 

(70) 

where 

MFEEDJ  =  feedback  torque  about  the  y  body  frame  axis 
MFEEDK  =  feedback  torque  about  the  z  body  frame  axis 
The  gains,  K  and  C,  must  be  chosen  so  that  the  feedback 
torque  will  drive  the  angle  of  attack  to  an  acceptably  small 
value  in  a  short  period  of  time.  They  can  not  be  too  large, 
though,  or  else  the  gimbal  limits  of  the  missile's  stage  1 
nozzle  will  be  exceeded.  The  gimbal  limits  for  a  MM  III 
missile  happen  to  be  about  40°.  By  starting  with  Euler's 


equations  with  a  restoring  torque  term  added,  a  simplified  set 
of  equations  can  be  solved  to  produce  approximate  values  for 
the  gains.  Euler’s  equations  are: 


xxP 

+ 

1  yy 

lzz^z  =  0 

(71) 

yyq 

+ 

(  Ixx 

1zz>Pz  -  Fxze  -  0 

(72) 

ZZZ 

+ 

(Iyy 

^x^PS  -  Fx*e  =  0 

(73) 

Since  lyy  =  I zz  for  an  axisymmetric  body,  Eq  (71)  reduces  to 
p  =  constant.  Next,  the  assumption  is  made  that  the  velocity 
vector  is  a  constant  on  the  timescale  of  the  coning  motion  so 
that  Eqs  (72)  and  (73)  become  a  set  of  two  equations  with  two 
unknowns.  This  assumption  means  that: 

i  ..  .  •  ••  . 

=  q/  ey=  q»  0t=  r,  f^=  r 

Using  this,  along  with  Eqs  (64)  and  (65),  Eqs  (72)  and  (73) 
change  to: 


d  + 

yyOy 

<Jxx  - 

Iyy)p®« 

-  Fx(Kfiy 

+  c6y)  =o 

(74) 

S  4- 

yy°z 

(Iyy 

Ixx^ 

-  Fx'KS 

+  c  g>z)  =0 

(75) 

Putting  Eqs  (74)  and  (75)  into  matrix  form: 


Iyy//px 

0 


)/F 


yy 


)/F 


o 


(76) 


Substituting  the  following  for  6  )  i 

r*t 

yields : 

i  XVfx  -  KX  -  c  -pX(ivv-ixx)/fx 

.  =0  (77) 

pX(Iyy-Ixx)/Fx  IyyX  /Fx  -  KX  -  C 

Which  produces  the  following  fourth  order  equation: 

Iyy2X4/Fx2  -  2KIyyX3/Fx  +  (p2(Iyy-Ixx)2/Fx2 

-  2CIyy/Fx  +  K2)Xl  +  2KCX  +  c2  =  0  (78) 

By  substituting  values  for  the  known  constants  Ixx,  Iyy, 

Fx,  and  p  a  fourth  order  equation  can  be  set  up  in  K  and  C. 
Routhe's  stability  criterion  (15:185-188)  can  be  used  to  show 
that  in  order  for  the  system  to  be  stable,  K  <  0.  By  means  of 
tr ial-and-er ror  it  was  found  that  for  K  =  -0.9  and  C  =  -10  the 
roots  are  X,^  =  -0.3019+  j3.139,  X3ij=  -0.1538  ±  jl.599. 

These  roots  are  all  stable  with  the  longest  damping  period 
being  about  six  seconds.  This  result  is  satisfactory  because 
for  shorter  damping  periods  the  missile  goes  through  wild 
oscillations  to  achieve  the  desired  angle  of  attack  causing 
the  aerodynamic  forces  to  take  control.  Figures  5.4  and  5.5 
show  the  missile's  response  without  the  feedback  controller 
and  Figures  5.6  and  5.7  show  the  same  test  case  run  with  the 


feedback  controller  installed.  Comparing  Figures  5.4  and  5.6, 
notice  that  the  angles  of  attack  around  the  y  and  z  body  axes 
are  reduced  quicker  with  the  controller  installed  than  without 
it,  which  is  the  result  desired.  A  comparison  of  Figures  5.5 
and  5.7  shows  that  without  the  feedback  controller  the  total 
angle  of  attack  keeps  oscillating  even  after  a  long  period  of 
time  but  with  the  controller  installed,  the  angle  of  attack 
goes  to  zero  in  about  25  seconds.  A  decaying  exponential 
curve  is  superimposed  on  the  angle  of  attack  plot  in  Figure 
5.7.  The  exponential  curve  was  generated  by: 

<*  =  oioexp(-Xt)  (79) 


where 

ot  =  angle  of  attack 
c*0  =  initial  kick  angle 

X  =  the  longest  decay  constant  =  0.1538 
Notice  that  the  simulation  stays  well  below  the  predicted 
angle  of  attack  during  the  first  25  seconds.  This  is  probably 
due  to  the  jet  damping  terms  which  help  damp  out  the  angle  of 


attack . 


Flight  path  angle  (radians) 
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Chapter  Six  Results  and  Conclusions 


Position  -  Velocity  Errors 


The  first  noticeable  difference  in  the  comparison  of  the 
trajectories  of  a  spinning  and  nonspinning  missile  is  that  the 
stage  1  burnout  state  vectors  are  different.  Table  6.1  lists 
the  burnout  position  and  velocity,  relative  to  the  launch 
point,  for  three  cases,  namely,  15°  initial  kick  angle  with 
and  without  symmetry  axis  spin  and  16.5°  kick  angle  with 


spin . 


Table  6.1  Position  -  Velocity  Comparisons 


X,  ft 

Y,  ft 
Z  ,  ft 

u,  ft/sec 

v,  ft/sec 

w,  ft/sec 


15°,  no  spin 
131,680 


-80,106 

5,086 


-3,723 


|Vcm| ,  ft/sec  6,303 


15°,  spin 
134,650 
10,418 
-72,713 
5,246 


-3,394 

6,267 


16.5°,  spin 
131,760 
11,072 
-78,954 
5,092 


-3,673 

6,300 


In  the  last  two  cases  the  spin  rate  was  1.571  rad/sec.  The 
altitude,  represented  by  x  in  the  table,  differs  by  about  3000 
ft  between  the  spinning  and  nonspinning  cases  with  the  same 
initial  kick  angle.  By  spinning  the  missile,  the  altitude 
achieved  from  the  same  initial  kick  angle  is  greater  but  the 


downrange  distance  is  less.  The  differences  between  the  y 
direction  position  and  velocity  is  explained  by  the  fact  that 
the  missile's  flight  path  describes  a  coning  motion  in  the 
spinning  case  which  causes  some  dispersion  in  the  y  direction 
This  could  probably  be  eliminated  by  introducing  a  small 
initial  kick  angle  in  the  negative  y  direction.  The  Vcm 
row  in  the  table  represents  the  magnitude  of  the  velocity 
vector  for  each  case.  Notice  that  the  16.5°,  spinning  case 
and  the  15°/  nonspinning  case  are  nearly  the  same.  The  x  and 
z  position  components  compare  favorably  for  these  two  cases 
also,  so  it  seems  that  spinning  the  missile  at  the  relatively 
slow  rate  of  1.571  rad/sec  forces  the  trajectory  analyst  to 
choose  a  slightly  larger  initial  kick  angle  in  order  to 
achieve  the  same  burnout  vector  as  the  nonspinning  case. 
Figure  6.1  shows  the  relationship  between  stage  1  burnout 
altitude  and  initial  kick  angle  in  the  same  flight  regime  as 
that  of  Table  6.1.  The  altitude  that  the  15°,  nonspinning 
case  achieves  corresponds  to  the  16.5°,  spinning  case  just  as 
Table  6.1  shows.  The  conclusion  is  that  there  is  almost  no 
performance  penalty  due  to  spin. 


Gimbal  Angles  | ^ 

The  stage  1  engine’s  expansion  nozzle  on  most  missiles  is 
gimballed  so  that  minor  corrections  can  be  made  to  the  flight 
path  during  the  boost  phase.  As  stated  before,  the  MM  III 
gimbal  angle  limit  is  about  40°  so  this  is  a  hardware 
limitation  that  must  be  considered  in  this  study.  Figure  6.2 
shows  a  "worst  case"  plot  of  gimbal  angle  variation  with  time. 

The  feedback  controller  produces  the  largest  restoring  moments 
during  the  first  few  seconds  of  flight  when  the  angle  of 
attack  is  the  largest  due  to  the  initial  kick  angle.  A  high 
spin  rate,  20  rad/sec,  was  chosen  with  an  initial  kick  angle 
of  15°  to  see  if  the  spin  rate  produced  gimbal  angles  that 


were  too  large.  The  plot  shows  that  even  at  high  spin  rates 


the  gimbal  angle  only  changes  about  11°  at  most.  This  is 
well  within  the  hardware  limits  of  the  missile  system.  The 
rate  at  which  the  gimbal  angle  changes  may  be  of  concern 
though.  The  gimbal  angle  changes  about  11°  in  0.4  seconds. 
This  may  be  a  concern  if  the  gimbal  actuator  can  not  respond 
at  this  rate. 

Conclusions 

The  model  tested  and  discussed  in  this  report  was 
subjected  to  varying  spin  rates  and  compared  to  nonspinning 
cases.  It  was  necessary  to  design  a  feedback  controller  to 
minimize  the  angle  of  attack  and  with  this  modification 
included  it  seems  to  be  possible  to  spin  a  missile  without 
adverse  effects  to  the  trajectory.  Some  adjustments  in  the 
trajectory  must  be  made  to  produce  the  same  stage  1  burnout 
state  vector  as  the  nonspinning  model  but  the  variations  are 
predictable  so  the  burnout  vector  is  also  predictable. 


Appendix  A 


A  listing  o£  the  program  which  was  used  for  this  study  is 
shown  below.  It  was  run  with  the  input  file  shown  in  Appendix 
B  on  a  VAX  11/785  computer  with  a  VMS  operating  system. 


program  misspin 


i************* 


program  misspin  is  a  simple  dynamics  propagator 
it  takes  input  state  vector  y  and  propagates 
it  from  to  to  tf. 

written  by:  Capt  Robert  W.  Bandstra  AFIT/GA-85D 

yy(i)  is  an  array  to  read  the  a  matrix  into, 
t  is  the  time 
to  is  the  initial  time 

tf  is  the  final  time 

y(i,j)  is  the  array  that  contains  the  state  vectors 

n  is  the  number  of  equations 

h  is  the  timestep 

mode  is  not  currently  used 

nl  is  the  number  of  times  the  outside  loop  is  performed 

n2  is  the  number  of  times  the  inside  loop  is  performed 

nstp  is  used  as  an  input  to  vary  the  step  size 

common  /ham/  t  ,y ( 18 , 4)  ,  f ( 18 , 4)  ,er r ( 18)  ,n , h ,mode 
double  precision  t ,y , f ,err  ,h  ,tf 

common  /miscon/  1 1 , 1 2 , 1 3 , 14 ,m0 ,m2 ,m3 ,m4 , rad , r 2 ,r 3 ,r 4 , 
1  bdot ,mtop , xtop , xbot ,alf a , f lpa 
double  precision  1 1 , 1 2 , 1 3 , 1 4 , mO ,m2 , m3 ,m4 , r ad , r 2 , r 3 , r 4 
1  bdot , mtop ,xtop ,xbot ,alf a , f lpa 
common  /count/  j 

common  /feedbk/  mf eed j , mf eedk ,xe ,ye , ze , f x , f y , f z , 

1  theta j  , thetak 

double  precision  mf eed j ,m£ eedk ,xe ,ye , ze , f x , f y , f z , 

1  theta j  ,  thetak 
real*8  yy(18) 

j  =  0 

j  is  a  switch  that  makes  sure  a  section  of  atmos  is  only 
executed  once. 


open  input  file 


open ( 2 , f ile  = 1  input .dat2'  ,status  =  'old') 
open  output  files 


open ( 3  ,  f ile= ' output .dat2'  ,status='old') 
open (4,file='plot.dat3' ,status='old' ) 
open (7 , f ile= 1  plot .dat4'  ,status='old') 

read  input 

read  (2,10) to ,  tf 
read  (2,20) (y(i,l)  ,i  =  l,3) 
read  (2,20)  (y(i,l)  ,i  =  4,6) 
read  (2,20) (y ( i , 1) ,i=7,9) 
read  ( 2 , 20)  (yy ( i )  , i  =  10 , 12) 
read  (2 ,20) (yy ( i) ,i  =  13 , 15) 
read  ( 2 , 20)  (yy ( i ) , i  =  16 , 18) 

10  format (2x ,2e20 . 13) 

20  for mat ( 2x , 3e20 . 13 ) 

dir.  cos  must  be  converted  from  y  (angles)  to  cos(y) 

do  25  i  =  10,18 
y(i ,1)  =  dcosd (yy (i) ) 

25  continue 

y  (12,1)  =  -y (12,1) 

read  mode,  number  of  steps  for  numerical  integration, 
and  plot  flag 

read  ( 2 , 30) mode ,nstp , iplot 
30  f  or  mat  ( 2x  ,  i  5 ,  i  6  ,  i  5) 

print  ics 

50  format (//, 2x , 'dynamics  propagator ',//, 2x ,' to ,  tf:' 

1  2e20. 13 ,//, 2x , 1  initial  state  vector', 

2  /,2x,3e20.13,/,2x,3e20.13,/,2x,3e20.13,/, 

3  /,2x,3e20.13,/,2x,3e20.13,/,2x,3e20.13,/, 

4  2x,'mode,  stps : ' , 2i 6 ,/) 
write(3,50)to,tf,(y(i,l) ,i=l,18) , mode ,nstp 
call  strtup 

setup  rest  of  haming  initialization 

number  of  odes: 
n  =  18 
t  =  to 

timestep  calculation 

nl  sets  the  number  of  points  sent  to  the  output  files 


nl  =  minO ( 244 , nstp) 
n2  =  1  +  nstp/nl 
h  =  (tf  -  to)/(nl*n2) 
c 

c  initialize  haming 

c 

nxt  =  0 

c  haming  initialization  flag 
call  haming (nxt ) 
if(nxt.eq.O)  stop  99 
c 

c  numerical  integration  loop  -  one  timestep  per  call 

c 

do  1000  int  =  l,nl 
do  900  nit2  =  l,n2 
call  haming (nxt) 

900  continue 

wr ite ( 3 , 61)  t,alfa, flpa, thetaj, thetak 
61  format(2x,'t  = '  ,el6 . 9  ,/ , 2x  , 

2  ' alf a  ='  ,e20 . 13  , '  flpa  = '  ,e20 . 13  ,/ ,  2x  , 

3  'thetaj  =',el6.9,'  thetak  =',el6.9) 

c  wr ite ( 3 ,60) t ,  (y ( i  ,nxt )  ,i=l,18)  ,alfa , flpa , thetaj , thetak 

60  format (/, 2x ,' t  = ' ,e20 . 13 ,/ , 2x , ' y  = '  , 3e20 . 13  ,/ ,  2x  , 

1  3e20.13,/,2x,3e20.13,/,2x,3e20.13,/,2x,3e20.13,/,2x, 

2  3e20.13,/,2x,'alfa  =',e20.13,'  flpa  = '  ,e20 . 13  ,/ ,  2x  , 

3  'thetaj  *'fel6.9»'  thetak  =',el6.9) 
write(4,95)  t,alfa 

write(7f95)  t,flpa 
95  f or  mat (2x , 2el0 . 3) 

if ( iplot .ge . 0)  write(4,90)  y ( 8  ,  nxt )  ,y ( 9  ,nxt ) 

90  f or  mat ( 2x , 2e20 . 13 ) 

1000  continue 
c 

c  print  final  position 

c 

write(3,70)  (y(i,nxt)  ,i  =  l,18) 

70  f or  mat (/ , 2x , ' at  time  tf  2x  state  vector  is',/, 

1  2x,e20.13,/,2x,e20.13,/,2x,e20.13,/ 

2  2x,e20.13,/,2x,e20.13,/,2x,e20.13,/ 

3  2x,e20.13,/,2x,e20.13,/,2x,e20.13,/ 

4  2x,e20.13,/,2x,e20.13,/,2x,e20.13,/ 

5  2x,e20.13,/,2x,e20.13,/,2x,e20.13,/ 

6  ,2x,e20.13,/,2x,e20.13,/,2x,e20.13,/) 
close (unit=2) 

close (unit=3) 
close (unit=4) 
close (unit=7) 
stop 
end 


c 

subroutine  strtup 
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c 

c* 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 


set  up  missile  parameters  at  beginning  of 


flight 


tmass  is  the  total  mass  of  the  missile  at  launch  in  slugs 
rad  is  the  radius  cf  stage  i  in  ft. 
r2  is  the  radius  of  stage  ii  in  ft. 

r3  is  the  radius  of  stage  iii  in  ft. 

r4  is  the  radius  of  stage  iv  in  ft. 

11  is  the  length  of  stage  i  in  ft. 

12  is  the  length  of  stage  ii  in  ft. 

13  is  the  length  of  stage  iii  in  ft. 

14  is  the  length  of  stage  iv  in  ft. 

tlen  is  the  total  length  of  missile  at  launch  in  ft. 
gm  is  the  grav.  parameter  of  the  earth  in  ft**3/sec**2 
thrust  is  the  thrust  of  stage  i  in  lbf. 
mO  is  the  initial  mass  of  stage  i  in  slugs 

m2  is  the  mass  of  stage  ii  in  slugs 

m3  is  the  mass  of  stage  iii  in  slugs 

m4  is  the  mass  of  stage  iv  in  slugs 

mtop  is  the  sum  of  m2,  m3,  &  m4  in  slugs 

xtop  is  the  location  of  the  com  of  the  top  3  stages  in  ft 
maspro  is  the  init.  mass  of  the  stage  1  prop,  in  slugs 
mdot  is  the  mass  flow  rate  of  stage  i  in  slug/sec 
burnt  is  the  stage  1  engine  burn  time  in  sec. 
bdot  is  the  rate  of  change  of  inner  radius 
of  stage  1  propellant  in  ft/sec 
bf  is  the  final  inner  rad.  of  stage  1  at  end  of  burn  in  ft 
rho2sq  is  the  square  of  the  offset  distance  of  the  center 
of  mass  flow  in  ft**2 


common  /misdat/  gm, thrust , mdot , tmass , burnt , maspro , mass , 

1  xxi,yyi,zzi,dixxdt ,d iyydt ,dizzdt,rho2sq,xbar,gainl, 

2  gain2 ,ga in3 ,gain4 

double  precision  gm, thrust , mdot , tmass , burnt , maspro , 

1  mass,xxi,yyi,zzi,dixxdt ,d iyydt ,d i zzdt , rho2sq  ,xbar  , 

2  gainl ,gain2 ,gain3  ,gain4 

common  /miscon/  1 1 , 1 2 , 1 3 , 1 4 , mO , m2 ,m3 ,m4 , r ad , r 2 , r 3 , r 4 , 

1  bdot , mtop , xtop ,xbot ,alf a , f lpa 
double  precision  11 , 1 2 , 13 , 1 4 , mO , m2 , m3 , m4 , r ad , r 2 , r 3 , r 4 , 

1  bdot ,mtop , xtop , xbot , alf a , f lpa 
double  precision  tlen ,bf  ,x2  , x3  ,x4 

read  missile  parameters 


read (2,10) 
read (2,10) 
read (2,10) 
read (2,10) 
read (2 ,10) 
read (2,10) 
read (2,11) 


tmass, rad, tlen 

thrust , mdot ,burnt 

11,12,13 

14 ,m0 ,m2 

m3 ,m4 , maspro 

r2,r3,r4 

gainl  ,gain2 


a 
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read (2,11)  gain3,gain4 
0  format(2x,3e20.13) 

1  f or  mat ( 2x , 2e20 . 13 ) 

print  missile  parameters 

write(3,20)  tmass,r ad, tlen, thrust ,mdot , burnt, 11, 12, 13, 

1  14,m0,m2,m3,m4,maspro,r2,r3,r4,gainl,gain2,gain3,gain4 
0  f or  mat (/ , 2x , ' tmass  =',el6.9,'  rad  =  1  ,el6. 9 ,/ ,  2x , 

l'tlen  =',el6.9,'  thrust  = ' , el6 . 9 ,/ , 2x , ' mdot  =',el6.9, 

2'  burnt  = •  ,el6 . 9 ,/ , 2x , ' 11  =*,el6.9,'  12  = 1  ,el6 . 9 ,/ , 2x , 

3 ' 1 3  = ' ,el6 . 9 , '  14  = ' ,el6 . 9 ,/ , 2x , ' mO  =',el6.9,'  m2  =’, 

4  el6 • 9 ,/ , 2x , '  m3  =',el6.9,'  m4  = ' ,el6.9 ,/ , 2x , ' maspro  =', 

5  el6 • 9 , '  r 2  = • ,e 16 . 9 ,/ , 2x , ' r 3  =',el6.9,'  r4  =',el6.9, 

6  /,2x,'gainl  =',el6.9,'  gain2  =',el6.9, 

7  /,2x,'gain3  =',el6.9,'  gain4  =',el6.9) 

gm  =  1.40764688 2d +16 

bf  =  rad  -  0.08d+00 
bdot  =  bf/burnt 
mtop  =  m2  +  m3  +  m4 
xbot  =  ll/2.0d+00 
x2  =  11  +  1 2/2 . 0d+00 
x3  =  11  +  12  +  1 3/2 . 0d+00 
x4  =  11  +  12  +  13  +  14/2. 0d+00 
xtop  =  (x2*m2+x3*m3+x4*m4)/mtop 
rho2sq  =  rad**2 .d+00/2 .d+00 
write(3,30)  bf ,bdot , mtop , xbot ,xtop ,rho2sq 
0  format ( 2x bf  =',e20.13,'  bdot  = ' ,e20 . 13 ,/ , 2x , ' mtop  =' 

1  ,e20.13,'  xbot  = '  ,e20 . 13  ,/ , 2x , ’ xtop  =',e20.13, 

2  '  rho2sq  =',e20.13) 
return 

d 


subroutine  haming(nxt) 


haming  is  an  ordinary  differential  equations  integrator 
it  is  a  fourth  order  predictor-corrector  algorithm 
which  means  that  it  carries  along  the  last  four 
values  of  the  state  vector ,  and  extrapolates  these 
values  to  obtain  the  next  value  (the  prediction  part) 
and  then  corrects  the  extrapolated  value  to  find  a 
new  value  for  the  state  vector. 


the  value  nxt  in  the  call  specifies  which  of  the  4 
values  of  the  state  vector  is  the  'next*  one. 
nxt  is  updated  by  haming  automatically,  and  is  zero  on 
the  first  call 


the  user  supplies  an  external  routine  rhs(nxt)  which 


■m 


evaluates  the  equations  of  motion 

common  /ham/  x  ,y  ( 18 , 4)  ,  f ( 18 , 4)  ,er rest ( 18)  ,n ,h , mode 
double  precision  x  ,y ,  f ,er res t  ,h  ,hh  ,xo 

all  of  the  good  stuff  is  in  this  common  block, 
x  is  the  independent  variable  (  time  ) 
y(18,4)  is  the  state  vector-  4  copies  of  it,  with  nxt 
pointing  at  the  next  one 

f  (18,4)  are  the  equations  of  motion,  again  four  copies 
a  call  to  rhs(nxt)  updates  an  entry  in  f 
errest  is  an  estimate  of  the  truncation  error 
-  normally  not  used 

n  is  the  number  of  equations  being  integrated 
h  is  the  time  step 

mode  is  0  for  just  eom,  1  for  both  eom  and  eov 
to 1  =  0.0000000001 

switch  on  starting  algorithm  or  normal  propagation 
if  (nxt)  190,10,200 

this  is  hamings  starting  algorithm. .. .a  predictor  - 
corrector  needs  4  values  of  the  state  vector,  and  you 
only  have  one  -  the  initial  conditions, 
haming  uses  a  picard  iteration  (slow  and  painfull)  to 
get  the  other  three. 

if  it  fails,  nxt  will  still  be  zero  upon  exit, 
otherwise  nxt  will  be  1,  and  you  are  all  set  to  go 

10  xo  =  x 

hh  =  h/2 . 0d+00 
call  rhs(l) 
do  40  1  =  2,4 
x  =  x  +  hh 
do  20  i  =  l,n 

20  y  ( i  ,1 )  =  y  ( i  ,  1-1)  +  hh*f(i,l-l) 
call  rhs(l) 
x  =  x  +  hh 
do  30  i  =  l,n 

30  y ( i  ,  1)  =  y (i ,1-1)  +  h*f(i,l) 

40  call  rhs  ( 1 ) 
jsw  =  -10 
50  isw  =  1 

do  120  i  =  l,n 

hh  =  y ( i , 1)  +  h*  (  9 . 0d  +  00* f ( i , 1)  +  19 . 0d  +  00* f ( i  ,  2) 


+  h*f (i ,1) 


hh  =  y ( i , 1)  +  h*  (  9 . 0d  +  00* f ( i  ,  1)  +  19 . 0d  +  00*f ( i  ,  2) 

1  -  5 • 0d+00*  f ( i , 3 )  +  f ( i  ,  4)  )  /  24.0d  +  00 

if(  dabs (  hh  -  y ( i , 2 )  )  .It.  tol  )  go  to  70 
isw  =  0 

70  y (i ,2)  =  hh 

hh  =  y ( i , 1)  +  h* (f  (i  ,1) +4.0d  +  00*f (i  ,2) +f (i ,3) )/3.0d  +  00 
if  (dabs (hh-y ( i , 3) )  .It.  tol)  go  to  90 
isw  =  0 

90  y ( i , 3 )  *  hh 


•  •  *1 
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hh  =  y  (  i , 1)  +  h*(  3  .  Od  +  OO*  £  (  i  ,  1)  +  9  .  Od  +  OO*  f  ( i  ,  2 ) 

1  +  9  .  Od  +  OO*  f  (  i  ,  3)  +  3  .  Od  +  OO*  f  (  i  ,  4 )  )  /  8. Od  +  OO 

if (  dabs (hh-y ( i , 4) )  .It.  tol  )  go  to  110 
isw  =  0 

110  y ( i  , 4)  =  hh 
120  continue 
x  =  xo 

do  130  1  =  2,4 
x  =  x  +  h 
130  call  rhs(l) 

if(isw)  140,140,150 
140  jsw  =  jsw  +  1 

if (jsw)  50 ,280 , 280 
150  x  =  xo 
isw  =  1 

jsw  =  1 

do  160  i  =  l,n 
160  errest(i)  =  0.0 
nxt  =  1 

go  to  280 

190  jsw  =  2 

nxt  =  iabs(nxt) 

c 

c  this  is  hamings  normal  propagation  loop  - 

c 

200  x  =  x  +  h 

npl  =  mod  (nxt,  4)  +  1 
go  to  ( 210 ,230)  ,  isw 
c  permute  the  index  nxt  modulo  4 

210  go  to  (270,270,270,220) , nxt 
220  isw  =  2 

230  nm2  =  mod (npl, 4)  +  1 

nml  =  mod(nm2,4)  +  1 

npo  =  mod (nml, 4)  +  1 

c 

c  this  is  the  predictor  part 

c 

do  240  i  =  l,n 

f(i,nm2)  =  y(i,npl)  +  4.0d+00*h*(  2 . 0d+00* f ( i ,npo) 

1  -  f ( i ,nml)  +  2 . 0d+00*f ( i ,nm2)  )  /  3. 0d+00 

240  y ( i ,npl)  =  f(i,nm2)  -  0 . 9 25619835*e r rest  (  i ) 
c 

c  now  the  corrector  -  fix  up  the  extrapolated  state 

c  based  on  the  better  value  of  the  equations  of  motion 
c 

call  rhs(npl) 
do  250  i  =  l,n 

y(i,npl)  =  ( 9. 0d+00*y ( i ,npo)  -  y(i,nm2)  +  3.0d+00*h* 

1  ( f ( i  ,npl)  +  2.0d  +  00*f (i,npo)  -  f (i , nml) ) )/8 . 0d  +  00 

errest(i)  =  f(i,nm2)  -  y(i,npl) 

250  y(i,npl)  =  y(i,npl)  +  0.0743801653  *  errest(i) 
go  to  ( 260 , 270)  , jsw 
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subroutine  rhs(nxt) 


c 

c*  *  * 


rhs  calculates  the  equations  of  motion  of  a  spinning 
missile  in  flight. 

y(l-3,nxt)  are  the  x,y,z  comp,  of  the  position  vector 

y(4-6,nxt)  are  the  x,y,z  comp,  of  the  trans.  vel.  vector 

y(7-9,nxt)  are  the  x,y,z  comp,  of  the  rot.  vel.  vector 

rcm  is  the  mag.  of  the  COM  position  vector  in  ft 

y(10-18,nxt)  are  the  dir.  cos  relationships  between  the 
body  frame  and  the  inertial,  earth-centered  frame, 
thrx  is  the  inertial  frame  x-dir.  thrust  in  lbf. 

thry  is  the  inertial  frame  y-dir.  thrust  in  lbf. 

thrz  is  the  inertial  frame  z-dir.  thrust  in  lbf. 

common  /ham/  t  ,y ( 18 , 4 )  ,  f ( 18 , 4)  ,er rest ( 18)  , n , h , mode 

double  precision  t  ,y  ,  f  ,er rest  ,h 

common  /misdat/  gm, thrust ,mdot ,tmass , burnt ,maspro , mass 

1  xxi,yyi,zzi,dixxdt  ,diyydt  ,d izzdt ,  rho2sq  ,xbar ,gainl , 

2  gain2  ,gain3 ,gain4 

double  precision  gm, thrust ,mdot ,tmass , burnt ,maspro , 

1  mass,xxi,yyi,zzi,dixxdt ,diyydt ,dizzdt,rho2sq,xbar , 

2  ga ini  ,gain2 ,gain3 ,gain4 
common  /aero/  xa ,ya , za , 1 ,m ,n 1 
double  precision  xa ,ya , za , 1 , m, nl 

common  /feedbk/  mf eed j , mf eedk , xe ,ye , ze , f x , f y , f z , 

1  theta j , thetak 

double  precision  mf eed j , mf eedk ,xe ,ye , ze , fx , fy , fz , 

1  theta j  ,  thetak 

double  precision  rcm, thrx , thry , thr z 
call  mominr (nxt) 
call  feedbk(nxt) 

thrx  =  y(10,nxt)*fx  +  y(13,nxt)*fy  +  y(16,nxt)*fz 

thry  =  y(ll,nxt)*fx  +  y(14,nxt)*fy  +  y(17,nxt)*fz 

thrz  =  y(12,nxt)*fx  +  y(15,nxt)*fy  +  y(18,nxt)*fz 

rcm  =  dsqrt (y (1 , nxt) **2 .d+00  +  y (2  ,nxt) **2  .d+00 
1  +  y ( 3 ,nxt ) ** 2 .d+00) 
f (1 ,nxt)  =  y (4  ,nxt) 
f (2  ,nxt)  =  y ( 5, nxt) 
f ( 3  ,nxt)  =  y ( 6  ,nxt ) 

f(4,nxt)  =  -gm*y ( 1 ,nxt ) /rcm** 3 . 0d+00  +  thrx/mass 
1  +  xa/mass 

f(5,nxt)  =  -gm*y ( 2 ,nxt ) /rcm** 3 . 0d+00  +  thry/mass 
1  +  ya/mass 

f(6,nxt)  =  -gm*y (3 ,nxt) /rcm**3 . 0d+00  +  thrz/mass 


1  +  za/mass 

£(7,nxt)  =  (  (yy i-zzi ) /xxi ) *y ( 8 ,nxt ) *y ( 9 ,nxt ) 

1  +  (mdot*xe* (ye*y (8 ,nxt) +ze*y (9 ,nxt) ) ) /xxi 

2  -  y (7  ,nxt) *dixxdt/xxi 

3  +  1/xxi  -  mdot*rho2sq*y (7 ,nxt) /xxi 
f(8,nxt)  =  (  (zz i-xxi ) /yy i ) *y ( 7 ,nxt ) *y ( 9 ,nxt ) 

1  -  (mdot*y ( 8 ,nxt) *xe**2 . Od+OO) /yyi 

2  -  y ( 8  ,nxt ) *d iyydt/yy i 

3  +  m/yyi 

4  +  mfeedj/yyi 

f(9,nxt)  =  ( (xxi-yy i ) /zzi) *y ( 7 ,nxt ) *y ( 8 ,nxt ) 

1  -  (mdot*y (9 ,nxt) *xe**2 . Od+OO) /zzi 

2  -  y ( 9  ,nxt ) *di zzdt/zzi 

3  +  nl/zzi 

4  +  mfeedk/zzi 

£(10,nxt)  =  y  (9  ,nxt) *y (13 ,nxt)  -  y ( 8  ,nxt ) *y ( 16 ,nxt) 
f(ll,nxt)  =  y (9  ,nxt) *y (14 ,nxt)  -  y (8 ,nxt) *y (17  ,nxt) 
f(12,nxt)  =  y ( 9  ,nxt ) *y ( 15 ,nxt )  -  y ( 8  ,nxt ) *y ( 18  ,nxt ) 
£(13,nxt)  =  -y (9 ,nxt) *y ( 10 ,nxt)  +  y ( 7  ,nxt) *y ( 16  ,nxt ) 
f(14,nxt)  =  -y ( 9  ,nxt ) *y ( 11 ,nxt )  +  y  (7  ,nxt) *y (17 ,nxt) 
f(15,nxt)  =  -y (9 ,nxt) *y (12 ,nxt)  +  y  ( 7  ,nxt) *y ( 18 ,nxt ) 
f(16,nxt)  =  y  (8,nxt) *y (10,nxt)  -  y ( 7  ,nxt ) *y ( 13 ,nxt ) 
f(17,nxt)  =  y (8 ,nxt) *y (11 ,nxt)  -  y ( 7  ,nxt) *y (14 ,nxt ) 
f(18,nxt)  =  y  (8,nxt) *y (12,nxt)  -  y ( 7 ,nxt ) *y ( 15 ,nxt ) 
return 


subroutine  mominr(nxt) 


mominr  computes  the  moments  of  inertia  and  their  first 
derivatives  as  a  function  of  time. 

xxi  =  mom.  of  iner.  about  the  body  x  axis  in  slug-ft**2 
yyi  =  mom.  of  iner.  about  the  body  y  axis  in  slug-ft**2 
zzi  =  mom.  of  iner.  about  the  body  z  axis  in  slug-ft**2 
dixxdt  is  the  der.  of  xxi  w.r.t.  time  in  slug-f t**2/sec 
diyydt  is  the  der.  of  yyi  w.r.t.  time  in  slug-f t**2/sec 
dizzdt  is  the  der.  of  zzi  w.r.t.  time  in  slug-f t**2/sec 
ml  is  the  instantaneous  stage  i  mass  in  slugs 
xbar  is  the  distance  between  the  com  and  the  bottom  of 
the  missile  in  ft. 

dl  is  the  moment  arm  distance  for  stage  i  in  ft. 

d2  =  moment  arm  dist.  for  the  rest  of  the  missile  in  ft 

cdO  is  the  base  drag  coefficient  (nondim) 

cdalfa  is  the  drag  coeff.  variance  with  angle  of  attack 

vcm  is  the  speed  of  the  com  in  ft/sec 

flpa  is  the  flight  path  angle  of  the  missile  in  rad. 

alfa  is  the  angle  of  attack  in  rad. 

cd  is  the  drag  coefficient  (nondim) 

pi  is  not  something  you  eat 


c  af  is  the  frontal  area  in  ft**2 

c  as  is  the  missile  surface  area  in  ft*2 

c  area  is  the  effective  aerodynamic  area  in  ft**2 

c  rho  is  the  air  density  in  slug/ft**3 

c  aaa  is  a  convenient  grouping  of  terms  used  later 

c  xcp  is  the  distance  of  the  center  of  pressure  from  the 

c  base  of  the  missile  in  ft. 

c  d  is  the  moment  arm  for  the  aero,  moment  in  ft 

c  xa  is  the  aero,  force  in  the  inertial  frame  x-dir  in  lbf 

c  ya  is  the  aero,  force  in  the  inertial  frame  y-dir  in  lbf 

c  za  is  the  aero,  force  in  the  inertial  frame  z-dir  in  lbf 

c  1  is  the  aero  moment  around  the  body  fr.  x-axis  in  lbf-ft 

c  m  is  the  aero  moment  around  the  body  fr.  y-axis  in  lbf-ft 

c  nl  =  aero,  moment  around  the  body  fr.  z-axis  in  lbf-ft 

c 

common  /ham/  t  ,y ( 18 , 4)  , f (18 , 4 ) ,er rest ( 18) ,n ,h , mode 
double  precision  t ,y  ,  f  ,er rest  ,h 

common  /misdat/  gm, thrust , mdot ,tmass , burnt ,maspro , mass , 

1  xxi ,yyi ,zzi ,dixxdt ,diyydt ,dizzdt ,rho2sq,xbar ,gainl , 

2  gain2 ,gain3 ,gain4 

double  precision  gm, thrust , mdot ,tmass , burnt ,maspro, 

1  mass , xxi ,yyi ,zzi ,dixxdt ,diyydt ,dizzdt ,rho2sq ,xbar  , 

2  gainl ,gain2  ,gain3 ,gain4 

common  /miscon/  1 1 , 1 2 , 1 3 ,1 4 ,m0 ,m2 ,m3 ,m4 , rad , r 2 ,r 3 ,r 4 , 

1  bdot,mtop,xtop,xbot,alfa,flpa 
double  precision  11 ,12 , 1 3 , 14 ,m0 ,m2 ,m3 ,m4 ,rad , r 2 ,r 3 , r 4 , 

1  bdot,mtop,xtop,xbot,alfa,flpa 
common  /aero/  xa ,ya , za , 1 , m,nl 
double  precision  xa ,ya,za ,1 ,m,nl 
double  precision  ml ,dl ,d2 ,cdO ,cdalf a , vcm,cd , 

1  pi,af,as,area,rho,aaa,xcp,d,alj 
mass  =  tmass-mdot*t 
if (t .gt .burnt)  mass  =  tmass-maspro 
ml  =  mO  -  mdot*t 

xbar  =  (xbot*ml  +  xtop*mtop) / (ml  +  mtop) 
dl  =  xbar  -  xbot 
d2  =  xtop  -  xbar 

xxi  =  (ml* (rad**2 .d+00-bdot**2 .d+00*t**2  .d  +  00)  + 

1  m2*r2**2.d+00  +  m3*r 3**2 .d+00  +  m4*r 4**2 .d+00) /2 .d+00 
dixxdt  =  (-2 .d+00*m0*t*bdot** 2 .d+00  -  mdot*rad**2 .d+00 

1  +  3.d+00*mdot*bdot**2.d+00*t**2.d+00)/2.d+00 
yyi  =  (ml* (3 .d+00*rad**2 .d+00  +  11**2. d+00)  +  m2* 

1  (3.d+00*r2**2.d+00+12**2.d+00)  +  m3* ( 3 .d+00*r 3**2 .d+00 

2  +  13**2. d+00)  +  m4* ( 3 .d+00*r 4** 2 .d+00  +  14**2. d+00)  - 

3  3 .d+00*ml*bdot** 2 . d+00* t** 2. d+00) /12. d  +  00 

4  +  mtop*d2** 2 .d+00  +  ml*dl** 2 .d+00 
diyydt  =  ( 2 .d+00*mdot/mass) * (mtop*d2* (-dl) 

1  +  ml*dl** 2 .d+00)  -  ( mdot/12 .d+00) * ( 3 .d  +  00*rad**2 .d  +  00 

2  +  11**2. d+00)  -  bdot**2 .d+00*m0*t/2 .d+00 

3  +  . 7 5d+00*bdot** 2 .d+00*mdot*t**2 .d+00-mdot*dl** 2 .d+00 
zzi  =  yyi 

dizzdt  =  diyydt 
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c  next,  calculate  the  aerodynamic  forces  and  moments 
c 

cdO  =  0 . 3d+00 
cdalfa  =  0.1d+00 

vcm  =  dsqrt (y ( 4 ,nxt) **2 .d+00  +  y ( 5 ,nxt) **2 .d+00  + 

1  y ( 6  ,nxt) **2 .d+00) 
flpa  =  dacos (y ( 4 ,nxt) /vcm) 

alj  =  dsqrt (y ( 10 ,nxt) **2 .d+00  +  y  (11  ,nxt) **2 .d+00  + 

1  y (12 ,nxt) **2 .d+00) 

alfa  =  dacos ( (y (4  ,nxt ) *y ( 10 ,nxt)  +  y ( 5 ,nxt) *y (11 ,nxt) 

1  +  y ( 6 ,nxt) *y (12 ,nxt) ) / (vcm*al j ) ) 
cd  =  cdO  +  cdalfa*alfa 
pi  =  3. 141592653589 8d +00 
af  =  pi*rad**2 .d+00 

as  =  pi* ( 2 .d+00* (rad*ll  +  r2*12  +  r3*13  +  r4*14) 

1  +  rad**2 .d+00) 

area  =  af *dcos (alf a)  +  as*dsin (alf a) 

zin  =  sngl (y (1 ,nxt) )  -  2 . 0925673e+07 

call  atmos (zin , dens) 

rho  =  dble(dens) 

aaa  =  -0 . 5d+00*cd*rho*area* vcm 

xa  =  aaa*y(4,nxt) 

ya  =  aaa*y(5,nxt) 

za  =  aaa*y(6,nxt) 

xcp  =  28.704d+00 

d  =  xbar  -  xcp 

1  =  0.0d+00 

m  =  d*aaa*  (y (16,nxt) *y(4,nxt)  +  y ( 17 ,nxt ) *y ( 5 ,nxt )  + 

1  y (18,nxt) *y(6,nxt) ) 

nl  =  -d*aaa* (y (13  ,nxt ) *y ( 4 ,nxt)  +  y (14,nxt) *y (5,nxt)  + 
1  y (15,nxt) *y (6,nxt) ) 
return 
end 

c* ********************************************************** 
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subroutine  atmos (zin ,dout) 


this  subroutine  is  an  atmospheric  model  that  calculates 
pressure,  density,  and  temp  at  altitudes  up  to  300km. 
input  is  altitude,  zin(meters),  in  call  statement, 
subroutine  parameters  are  gas  const,  r ( joules/kg-deg-k) 
and  the  S.L.  values  of  the  grav.  const,  g 0 (m/sec-sq) , 
atmospheric  pressure  at  sea  level,  p0(n/m-sq), 
and  the  density,  d0(kg/m-cu). 

the  first  order  gravitational  constant  is  b(l/m). 
the  following  subscripted  variables  are  used: 
z(i)  is  the  altitude  (km) 

t(i)  is  the  molecular  temperature  (deg-k) 

p(i)  is  the  pressure  (n/m-sq) 

d(i)  is  the  density  (kg/m-cu) 


l(i)  is  the  thermal  lapse  rate  (deg-k/km) 
common  /count/  j 

dimension  z(22) ,t(22) ,p(22) ,d(22) ,1(22) 
if (3 .eq. 1)  goto  200 

j  =  1 

r  =  287.0 

gO  =  9.806 

p0  =  1.01325e+05 

dO  =  1.225 

b  =  3 . 139e-07 

units  change  form  ft  to  km 
zin  =  zin/3 . 281e-t-03 
data  (z(i)  ,  i*l ,  22) / 

1  00.00,  11.019,  20.063,  32.162,  47.350, 

2  52.43,  61.590,  80.00,  90.00,  100.00, 

3  110.00,  120.00,  150.00,  160.00,  170.00, 

4  190.00,  230.00,  300.00,  400.00,  500.00, 

5  600.00,  700.00/ 

data  ( t ( i ) ,i=l,22)/ 

1  288.10,  216.65,  216.65,  228.65,  270.65, 

2  270.65,  252.65,  180.65,  180.65,  210.65, 

3  260.65,  360.65,  960.65,  1110.65,  1210.65, 

4  1350.65,  1550.65,  1830.65,  2160.65,  2420.65, 

5  2590.65,  2700.65/ 

data  (p(i) ,i=l,22)/ 

1  1.000e+00,2.284e-01,5-462e-02,8.567e-03,1.095e-03 

2  5.823e-04,1.797e-04,1.024e-05,1.622e-06,2.980e-07 

3  7.220e-08,2.488e-08,5.000e-09,3.640e-09,2.756e-09 

4  1.660e-09,6.869e-10,1.860e-10,3.977e-ll,1.080e-ll 

5  3 . 400e-12  ,  1.176e-12/ 
do  100  i  =  1,22 

z(i)  =  z(i)  *  1000 
P ( i )  =  P ( i )  *  p0 
d (i)  =  p(i)/(r*t (i) ) 

100  continue 

do  200  i  =  1,21 

1 ( i )  =  (t (i+1) -t (i) )/(z (i+1) -z(i) ) 

200  continue 

do  300  i  =  1,21 

if (zin.ge .z (i+1) )  go  to  300 

i f (abs (1 ( i ) ) . It . 1 . 0e-05)  go  to  400 

ql  =  1+b*  ( (t(i)/l(i)  ) —  z (  i  )  ) 

q2  =  (ql*g0) / (r*l (i ) ) 

tout  =  t ( i ) +1 ( i ) * ( zin-z ( i ) ) 

q3  =  tout/t ( i ) 

q4  =  q3** (-q2 ) 

q5  =  exp(  (b*g0* (zin-z(i) ) )/(r*l (i)  )  ) 
q6  =  q4*q5 
pout  =  p(i)*q6 
q7  =  q2  +  1 

dout  =  d(i)*(q3**(-q7))*q5/515.4 
go  to  500 


'.’V 


400  tout  =  t  (  i  ) 

q8  =  (-1. ) *g0* (zin-z (i) ) * (l.-b/2. ) * (zin+z (i ) ) / (r *t (i ) ) 
pout  =  exp (q8) *p ( i ) 

c  units  change  from  kg/m**3  to  slug/ft**3  (density) 
dout  =  exp (q8) *d ( i ) /515 . 4 
go  to  500 
300  continue 
500  return 
end 

c 

subroutine  feedbk(nxt) 
c 

c 

c  feedbk  calculates  the  missile's  angle  of  attack  in  the 
c  body  frame  and  returns  a  restoring  moment  to  rhs 


c 

vi  =  missile 

vel.  in 

the 

body 

frame , 

x-dir.  in 

ft/sec 

c 

vj  =  missile 

vel.  in 

the 

body 

frame , 

y-dir.  in 

ft/sec 

c 

vk  =  missile 

vel.  in 

the 

body 

f  r  a  me  , 

z-dir.  in 

ft/sec 

c  thetaj  is  the  angle  of  attack  about  the  body  frame  y  axis 

c  thetak  is  the  angle  of  attack  about  the  body  frame  z  axis 

c  xe  is  the  offset  of  mass  flow  from  the  body  x  axis  in  ft 

c  ye  is  the  offset  of  mass  flow  from  the  body  y  axis  in  ft 

c  ze  is  the  offset  of  mass  flow  from  the  body  z  axis  in  ft 

c  fx  is  the  thrust  in  the  body  frame  x-dir.  in  lbf 

c  fy  is  the  thrust  in  the  body  frame  y-dir.  in  lbf 

c  fz  is  the  thrust  in  the  body  frame  z-dir.  in  lbf 

c  mfeedj  is  the  feedback  mom.  required  to  null  out  the  AOA 

c  about  the  body  fr.  y  axis 

c  mfeedk  is  the  feedback  mom.  required  to  null  out  the  AOA 
c  about  the  body  fr.  z  axis 

c  gainl  is  the  damping  constant  for  ye 

c  gain2  is  the  gain  coefficient  for  ye 

c  gain3  is  the  damping  constant  for  ze 

c  gain4  is  the  gain  coefficient  for  ze 

c 

common  /ham/  t  ,y ( 18 , 4)  , f ( 18 , 4)  ,er res t ( 18)  ,n ,h , mode 

double  precision  t  ,y  ,  f  ,er rest  ,h 

common  /misdat/  gm, thrust ,mdot ,tmass , burnt ,maspro, 

1  mass,xxi,yyi,zzi,dixxdt ,diyydt ,dizzdt ,  rho2sq  ,xbar  , 

2  gainl  ,gain2  ,gain3  ,gain4 

double  precision  gm, thrust ,mdot ,tmass , burnt ,maspro , 

1  mass,xxi,yyi,zzi,dixxdt,diyydt,dizzdt,rho2sq,xbar, 

2  gainl  ,gain2  , gain 3 ,gain4 

common  /feedbk/  mf eed j , mfeedk ,xe , ye , ze , fx , fy , fz , 

1  theta j  ,  thetak 

double  precision  mf eed j ,mf eedk ,xe ,ye , ze , f x , fy , f z , 

1  thetaj  , thetak 
*  double  precision  vi,vj,vk 
c 

xe  =  -xbar 
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vi  =  y (10 ,nxt) *y (4 ,nxt)  +  y ( 11 ,nxt ) *y ( 5 ,nxt ) 

1  +  y ( 12  ,nxt ) *y { 6 ,nxt ) 

vj  =  y  (13  ,nxt) *y ( 4 ,nxt)  +  y ( 14  ,nxt ) *y ( 5 ,nxt ) 

1  +  y  ( 15  ,nxt ) *y ( 6 ,nxt ) 

vk  =  y (16  ,nxt) *y ( 4  ,nxt)  +  y ( 17 ,nxt ) *y ( 5  ,nxt ) 

1  +  y (18  ,nxt) *y ( 6  ,nxt) 
thetaj  =  datan(vk/vi) 
thetak  =  datan(vj/vi) 
ye  =  gainl*y (9,nxt)  +  gain2*thetak 
ze  =  gain3*y (8 ,nxt)  +  gain4*thetaj 

if  (dabs (ye/xe) .ge . 0 . 8 391d+00) ye  =  dsign ( 17 . 619d+00 ,ye) 
if  (dabs ( ze/xe ) .ge . 0 . 8391d+00) ze  =  dsign ( 17 . 619d+00 , ze ) 
fy  =  thrust*ye/dsqrt (xe**2 .d+00+ye**2 .d+00+ze**2.d+00) 
fz  =  thrust*ze/dsqrt (xe**2 .d+00+ye**2 .d+00+ze**2 .d+00) 
fx  =  dsqrt (thrust**2 .d+00  -  fy**2.d+00  -  fz**2.d+00) 
mfeedj  =  ze*fx 
mfeedk  =  -ye*fx 
return 
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Appendix  B 


A  sample  of  the  input  file  used  to  support  program 
misspin  is  listed  below.  This  file  is  referred  to  in  line  42 
of  misspin  as  ' input. dat2'  and  contains  all  the  inputs 
required  to  run  the  program. 


to  tf 

0.0000000000000d+00  0. 610000000000 Od +02 

x  y  z 

2. 092572257000 Od +07  0 . 0000000000000d+00  0 . 0000000000000d+00 


u  v  w 

5. 000000000000 Od +01  0 . 0000000000000d+00  0 . 0000000000000d+00 


p  q  r 

1.570796370506 3d +00  0 . 0100000000000d+00  0 . 0100000000000d+00 

all  al2  al3 

0.1 50000000000 Od +02  0 . 9000000000000d+02  0 . 7 500000000000d+02 

a  21  a22  a  2  3 

0. 900000000000 Od +02  0 . 0000000000000d+00  0 . 9000000000000d+02 

a31  a32  a33 

0.7 50000000000 Od +02  0 .9000000000000d+02  0 . 1500000000000d+02 

mode  nstp  iplot 
+0000+15000-  1 

tmass  rad  tlen 

2. 362155800000 Od +03  2 . 7 500000000000d+00  5 . 9900000000000d +01 

thrust  mdot  burnt 

2. 020000000000 Od +05  2 . 3351510000000d+01  6 . 1000000000000d+01 

11  12  13 

2. 530000000000 Od +01  1 . 3100000000000d  +  01  7 . 1000000000000d  +  00 

14  mO  m2 

1. 440000000000 Od +01  1 . 5711444000000d+03  4 . 8237707000000d+02 

m3  m4  maspro 

2. 53 3 101 300000 Od +02  5. 5324175000000d+01  1 . 4244421000000d+03 


Appendix  L 


A  sample  of  the  output  file  produced  by  program  misspin 


is  listed  below.  This  file  is  referred  to  in  line  46  of 


misspin  as  1  output .dat2 '  and  contains  the  output  from  a 
standard  run  of  the  program.  A  large  number  of  iterations  has 
been  deleted  between  t  =  2.0  and  t  =  60.0. 


Dynamics  Propagator 

to,  tf:  0. 000000000000 Oe +00  0 . 6100000000000e+02 

Initial  State  Vector 

0. 209257225700 Oe +08  0 . 0000000000000e+00  0 . 0000000000000e+00 
0. 500000000000 Oe +02  0 . 0000000000000e+00  0 . 0000000000000e+00 
0. 157079637050 6e +  01  0 . 1000000000000e-01  0 . 1000000000000e-01 

0. 965925826289 le +  00  0 . 0000000000000e  +  00-0 . 258819045102 5e+00 
0. 000000000000 Oe +00  0 . 1000000000000e+01  0 . 0000000000000e+00 
0. 258819045102 5e +00  0 . 0000000000000e+00  0 . 9659258262891e+00 
mode,  stps:  0  15000 


tmass  =  0 . 236215580e+04  rad  =  0 . 275000000e+01 
tlen  =  0. 599000000e+02  thrust  =  0 . 202000000e+06 


mdot 


0. 23351510 Oe +02  burnt  =  0 . 610000000e+02 


11  =  0. 25300000 Oe +02  12  =  0 . 13 1 000000e+0 2 

13  =  0. 71000000 Oe +01  14  =  0 . 144000000e+02 

mO  =  0. 15711444 Oe +04  m2  =  0 . 482377070e+03 

m3  =  0. 253310130e+03  m4  =  0 . 553241750e+02 

maspro  =  0 . 142444210e+04  c2  =  0. 215000000e+01 

r  3  =  0. 21500000 Oe +01  r4  =  0 . 215000000e+01 

gainl  =-0 . 900000000e+00  gain2  =-0 . 100000000e+02 

gain3  =-0 . 900000000e+00  gain4  =-0 . 100000000e+02 

bf  =  0. 267000000000 Oe +01  bdot  =  0 . 4377049180328e-01 

mtop  =  0.79101137 5000 Oe +03  xbot  =  0 . 1265000000000e+02 

xtop  =  0. 365426523171 Oe +02  rho2sq  =  0 . 3781250000000e+01 

t  =  0. 25000000 Oe +00 

alf a  =  0. 118398410980 2e +00  flpa  =  0 . 113 5910610128e+00 
thetaj  «  0. 11065052 4e +00  thetak  =  0 . 424722783e-01 
t  =  0. 500000000e+00 

alfa  =  0. 134902716244 5e -01  flpa  =  0 . 1576384415400e+00 
thetaj  =  0. 13486906 le -01  thetak  =-0 . 301353252e-03 
t  =  0.7 5000000 Oe +00 

alfa  =  0.6806496272313e-01  flpa  =  0 . 1 63 32062057 55e+00 
thetaj  =-0. 94991864 8e -02  thetak  =-0 . 674028930e-01 
t  =  0 . 100000000e+01 

alfa  =  0. 10 808 37 27 138 le +00  flpa  =  0 . 1492277787722e  +  00 
thetaj  =  0.193839957e-01  thetak  =-0 . 1063 57830e+00 
t  =  0. 12500000 Oe +01 

alfa  =  0.1119415232809e+00  flpa  =  0 . 1282593230884e+00 
thetaj  =  0. 563485640e-01  thetak  =-0 . 969292474e-01 


t  =  0. 15000000 Oe +01 

alf a  =  0. 854521381570 4e -01  flpa  =  0 . 110 5389365479e+00 
thetaj  =  0. 65742803 6e -01  thetak  =-0 . 547463508e-01 
t  =  0 . 1 75000000e+01 

alf a  =  0.4292624643032e-01  flpa  =  0 . 1021869574191e+00 
thetaj  =  0. 40700527 Oe -01  thetak  =-0 . 136580155e-01 
t  =  0. 20000000 Oe +01 

alf a  =  0.111657 341299 9e -02  flpa  =  0 . 1040975924Q51e+00 
thetaj  =  0.827320623e-03  thetak  =  0 . 749851508e-03 
**  portion  of  output  deleted  for  the  sake  of  brevity  *** 
t  =  0. 60000000 Oe +02 

alfa  =  0. 11 8 34 7 07 27 66 5e -02  flpa  =  0 . 5760804791962e+00 
thetaj  =-0.487526051e-03  thetak  =-0 . 107838847e-02 
t  =  0. 60250000 Oe +02 

alfa  =  0. 162708479222 6e -02  flpa  =  0 . 5767944135272e+00 
thetaj  =-0. 14934812 3e -02  thetak  =  0 . 645693414e-03 
t  =  0. 60500000 Oe +02 

alfa  =  0.2458572367868e-02  flpa  =  0 . 5774826214918e+00 
thetaj  =-0. 17083999 5e -02  thetak  =  0 . 176803841e-02 
t  =  0.607 50000 Oe +02 

alfa  =  0.1542471668112e-02  flpa  =  0 . 5781562255459e+00 
thetaj  =-0. 816599337e-03  thetak  =  0. 130858163e-02 
t  =  0 . 610000000e+02 

alfa  =  0. 794214125723 Oe -03  flpa  =  0 . 5788422181625e+00 
thetaj  =  0. 72741351 8e -03  thetak  =-0 . 318819283e-0 3 


at  time  t f 


state  vector  is 
0. 210603675494 6e +08 
0. 104177288390 4e +05 
-0. 727133850594 8e +05 
0. 5246253808848e+04 
0. 487078046067 9e +03 
-0. 3393730020379e+04 
0.1516615939657e+01 
0.7507226603557e-02 
0. 980864939723 4e -02 
0. 836678517382 2e +00 
0.77 5936422123 7e -01 
-0. 5421699781792e+00 
-0. 3489327010687e+00 
0.8385295590313e+00 
-0. 418466424804 le +00 
0. 422155218684 6e +00 
0. 5393027028031e+00 


0.7 28654627 392 Oe +00 
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launch-centered,  nonrotating  earth,  inertial  reference  frame  and  the 
missile  body  reference  frame.  The  moments  of  inertia  and  aerodynamic 
forces  are  derived  and  presented.  A  feedback  controller  is  derived 
which  proved  to  be  a  necessary  addition  to  the  system  in  order  to 
reduce  the  angle  of  attack.  The  angle  of  attack  of  the  missile 
produced  adverse  effects  on  the  burnout  vector  without  the  feedback 
controller  but  the  effects  are  reduced  considerably  with  the 
controller  included.  Problem  areas  include  possibly  excessive  nozzle 
gimbal  rates  caused  by  the  feedback  controller  and  the  need  to  change 
the  initial  kick  angle  if  the  missile  is  spinning  in  order  to  achieve 
the  same  burnout  conditions  as  a  nonspinning  missile. 
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